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ABSTRACT 


This report was prepared for the Proceedings of the Fourteenth 
Stanstead Seminar, sponsored by the Department of Meteorology, 

McGill University, and the National Center For Atmospheric Research, 
and was written with the Stanstead audience in mind. The theme of 
the Seminar was "The interaction between objective analysis and 
initi alizati chi". We first review a duality between optimum interpolation 
and variational objective analysis. We use this duality to set up 
a variational approach to objective analysis which uses prior 
information concerning the atmospheric spectral energy distribution, 
in the variational problem. In the wind analysis example we study, 
the wind field is partitioned into divergent and nondivergent parts, 
and a control parameter governing the relative energy in the two 
parts is estimated from the observational data being analyzed by 
generalized cross validation, along with a bandwidth parameter. We 
then propose a variational approach to combining objective analysis 
and initialization in a single step. In a simple example of this 
approach, data, forecast, and prior information concerning atmospheric 
energy distribution is combined into a single variational problem. 

This problem has (at least) one bandwidth parameter, one partitioning 
parameter governing the relative energy in fast and slow modes, and 
one parameter governing the relative weight to be given to observational 
and forecast data. It may be possible to estimate good values of 
all three parameters from the observational and forecast data. 



1. INTRODUCTION 


In this work, we exploit a duality between optimum Interpolation and 
variational objective analysis to set up a certain variational approach to 
objective analysis, as well as to propose an approach to combining objective 
analysis and variational Initialization into a single step. 

In Section 2 we describe the duality. In Section 3 we describe the 

variational approach to objective analysis In the context of the estimation 

of vorticity and divergence from observed wind vectors. The duality of 
Section 2 Is used to suggest how data concerning the energy spectral 
distribution In the atmosphere, as studied by, e.g., Baer {1974, 1981), 

Kasahara (1976), Kasahara and Purl (1981), Stanford (1979), may be used to 
help choose the form of the variational problem to be solved. The 

variational problem we solve to estimate vorticity and divergence has one 
"bandwidth" parameter (related to the half power point of the equivalent low 
pass filter) and one "partitioning" parameter, representing the relative 
allocation of energy to the divergent and non-divergent part of the wind. 
Both can be estimated by generalized cross validation (GCV). 

In Section 4 we begin a synthesis of several Ideas •• 1) the duality of 
Section 2, 11) the Idea of partitioning the "signal" Into divergent and 

non-divergent parts (which generalizes to the Idea of partitioning the signal 
into slow and fast modes) and ill) the modified Kalman filter Ideas as 
proposed by Ghil and coworkers (1981). The result is an argument that (In 
principle) one can set up the problem of estimating initial conditions by 
combining i) the forecast, 1i) the observational data, ill) possibly certain 
physical constraints and iv) (partial) prior Information concerning 

atmospheric spectral energy distribution, into a single variational problem. 
This result is in some sense a converse of Phillips (1982b), who argues that 
normal mode initialization can be done as part of optimum interpolation. The 
resulting variational problem as we propose it will have (at least) one 
bandwidth parameter, one balance parameter controlling the relative weight to 
be given to forecast data and observational data, and one partitioning 
parameter, governing the relative energy in the "signal" assigned to fast and 
slow modes. We conjecture that these three (control) parameters can be 
estimated dynamically from the data by GCV. They can also be chosen by trial 
and error. Objective analysis and (linear) normal mode initialization are 
thereby combined in one step. The estimation of one (or more) balance 

parameter(s) from the data may possibly avoid the pitfalls due to 

misspecification of the Kalman filter statistics as noted by Phillips 

(1982a). 

Some of this work is joint with D. R. Johnson. 
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2. ON A DUALITY BETWEEN OPTIMUM INTERPOLATION AND VARIATIONAL OBJECTIVE 
ANALYSIS 


We Mill describe this duality for the analysis of a univariate variable on 
the sphere (say 500 mb height) although the result is completely general. 
Let P denote a point on the sphere end let h(P) be (say), the 500 mb height 
minus the global average 500 mb height at P. Let the observations 
yi» • • • Vn modelled as 

yi » h(Pi) + ci . (2.1) 


where the are supposed to be zero mean independent measurement errors with 
common variance = E 


We suppose E h(P) * 0 and h has a prior covariance b R(P,Q), 
E h(P)h(Q) = b R(P,Q). 

Then using standard results in multivariate analysis 
E(h(P)|yj...yn) = (R(P.Pi). ••«(?, ? n) )(Rn + nXI)-l 
= hx(P), say 


( 2 . 2 ) 


where X = o^/nb and R,, is the n x n matrix with ij th entry R(P-j,Pj). 
Considering the expression on the right of (2.2) as a function of P, which we 
denote by hx(P), it is easy to see that if X = 0, then ho(P) interpolates to 
the data exactly, ho(Pi) = yi . i = l,2,...n, whereas for X > 0, hx(P) smooths 
the data, and X controls the amount of smoothing -X is the "bandwidth" 
parameter, hx(P) evaluated at grid points is the "optimum interpolant" of 
Gandin, given that all the available data points are used simultaneously. 


Duality Theorem: Kimeldorf and Wahba (1970, 1971) Wahba (1978). 

For every covariance R(P,Q) satisfying //R^(P,Q)dPd0<“, there Is a 
variational problem for which hx(P) is the solution. It is: find h in Hr (a 

certain reproducing kernel Hilbert space) to minimize 


- I (yi - b(Pi) )^ + X J(h). 

n i=l (2.3) 

where J(h) is the square norm of h in Hr. 

We will now describe J(h) in - manner which we hope will make its 
meteorological usefulness clear. R(P,0) being a covariance, is a symmetric 
non-negative definite function and (given that it satisfies the hypothesis) 
has a so called Mercer-Hilbert-Schmidt expansion (Riesz -Sz.-Nagy (1955)) in 
<Hs eigenfunctions and eigenvalues 


•3* 




R(P.Q) * I ^ts Y£s(P) Yts(Q) . where J R(P,Q) Yis{Q)dQ • hs hs (P).{2‘4) 
t.s 

(This expansion Is a generalization of the factorization of a covariance 
matrix z as I * ror' where rr' = I and D Is diagonal). We have deliberately 
used a notation to suggest that the eigenfunctions of R are spherical 
harmonics, but that Is by no means necessary. In any case 


J(h) = I 



Is ^is 


where h^s = I h(P) Yiis(P)dP. 


(2.5) 


(If h were a vector and R a matrix, then we would have J(h) = h'R*^h). As an 
example, if the are the spherical harmonics and X£j = [£(A+l)]-2m £^en by 
using the fact that the spherical harmonics are the eigenfunctions of the 
Laplaclan, A Y^j = -£(i+l) Y^c, it Is not hard to show that J(h) = 

//(A^hj^ dP. 

m 

More generall>, if | Ov[(t) (t+1 )]'' |-2, then 

v=o 


m 

0(h) = /I I Ov (-A)''h|2 dP. 
v=o 


( 2 . 6 ) 


More details may be found in Wahba (1981a, 1981b), the use of Hough functions 
instead of spherical harmonics is briefly described in Wahba (1981b). If 

E h(P)h(Q) = b R(P,Q) , (2.7) 

then h has a Karhunen-Loeve expansion 

h(P) = b I h£s Y£c(P) where bhju = / h(P) Y^g (P) dP , (2.8) 


and the h^g are independent zero mean random variables with E h^g hfg' = 
bX£g, i,s = t',s', = 0 otherwise. The {X£g}, that is, the relative energy 

distribution by wave numbers may be estimated from historical data and/or 
obtained (roughly) from theory, see Baer (1974, 1981), Kasahara (1976), 

Kasahara and Puri (1981), Stanford (1979). 


3. ESTIMATION OF DIVERGENCE AND VORTICITY FROM THE OBSERVED WIND FIELD USING 
VECTOR SPLINES ON THE SPHERE. PARTITIONING OF THE DATA INTO DIVERGENT 
AND NON-DIVERGENT PARTS 

Given observed wind data (ui,v-j) at point P^ , i = 1,2,. . . n, we estimate 
the vorticity and divergence as follows. The stream function and velocity 
potential are expanded in spherical harmonics 

L £ L £ 

’ (P) = I, I au »Js(P) . t (P) = I I b,s Yis(P) . 

£=1 S=-£ £=1 s=-£ 


(3.1) 
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Then (for given 6,X), we find {a^s> l>£s> to minimize 
1 n 1 1 34 2 

- I(---(Pi)+ (Pi)-ui)^ (3.2) 

n i»l a 34 acos4^ 3X 

1 n 1 S'? 13* 

+ - i( (Pi) + (Pi)-Vi)^ 

n i=l acos4i 3X a 34 

1 

+ X[Ji(T) + 'J2(*)3 
5 


where 

Ji(t) I *^2(*) * ^ * (3»3) 

Z=1 


Given estimates of the ajtj and b£s* estimates of the wind, vorticity and 
divergence are obtained analytically for any desired P. The resulting 
estimated wind field is called a vector spline on the sphere. X is the 
bandwidth parameter, it controls the partitioning of the data vector 
(ui, . . . Up, vj, . . . Vp)' into a part due to "signal'' and a part due to 
noise. 6 can be viewed as a signal partitioning parmeter, it controls the 
partitioning of the "signal" part of the data into a divergent part and a 
non-divergent part. 

A Monte Carlo study was carried out to test the effectiveness of the 
method and to determine whether good estimates of X and 6 could be obtained 
from the data by GCV. The results are reported in Wahba (1982), and in 
preparation. The weights X£g(^) and X£s(^) were adapted from the data 
collected by Stanford (1979). Realistically scaled model 500 mb wind fields 
were generated via a Karhunen-Loeve expansion in stream function and velocity 
potential, the resulting "true" wind vectors at 114 North American weather 
stations computed and a random 2.5 m/sec rms measurement error added to each 
computed wind component. Good recovery of winds, vorticity and divergence in 
an area covered by the data grid and extending a small amount past it was 
obtained using the estimated X and 6. (The individual a£s» b£g are not 

recovered from only North American data). The results are sensitive to both 
X and 6. It can be seen from the experiments that a poor value of 6 causes 
obvious edge effects and that fixing 6 at values which oversuppress the 
divergent part of the wind tend to cause increased errors in the estimation 
of the non-divergent part. 


4. ON A SINGLE VARIATIONAL PROBLEM FOR MERGING DATA, FORECAST AND PRIOR 
KNOWLEDGE OF ATMOSPHERIC SPECTRAL ENERGY DISTRIBUTION 

We claim that a set of "moderately" reasonable assumptions concerning the 
forecast errors and measurement errors, along with prior knowledge concerning 
the spectral energy distribution in the atmosphere, leads to an optimal" 
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Initial state obtained as the solution of a single (large) variational 
problem. This variational problem will have the same number of "unknowns" as 
the degrees of freedom in the forecast model, thus raising questions of 
numerical feasibility. Vie brush aside computational questions for the 
present, adopting the point of view that it is worthwhile to examine an 
"optimal" variational problem, and then ask, how close can one come to 
computing a reasonable approximation to it. (Some numerical shortcuts appear 
in Bates and Wahba (1982)). One of the advantages of examining the 
variational form is that it is . fairly evident how to add side constraints 
based on the physics. 

A (simplified) example goes as follows. Consider a spectral model where 
the state of the atmosphere is expressed in terms of Hough functions. Using 
notation similar to Tribbia (1982), and considering a single level, let 



Hnax 




^ax 

I 

k=l 


Vk Bk® 


(4.1) 


(U,V,i>) is the "true" wind and geopotential field, and and H|(® are 

rotational (slow) and gravitational (fast) Hough functions. We suppose that 
the true "state" of the atmosphere is adequately described by the 

vector e = (X:Y) = (Xi, . . . Xm : Yi, . . . Ym ) and an analysis 
■ max « max 

consists of obtaining an updated estimate 6 of 6, given a forecast 
eF . (jF:f) . (X,F \„=’'l'' 


given data (ui,vi,4»i) representing observations on the wind and geopotential 
height at point P-j , and given the "prior knowledge" concerning atmospheric 
energy distribution obtained, e.g. from data like that collected by Kasahara 
and Puri (1981). This prior knowledge is of the form 
2 2 


EX = b, x.R 



(4.3) 


where we assume that the relative energies XjR and X|^G are known but that bj 
and b 2 are not. We will assume that all cross covariances E X^Xj, E X-jYj 
and E Y^Yj can be approximated by 0. Suppose further that 


E(Xj-XjF)(x,,_X,,F) = . E(Yj-YjF)(Yk-YkF) = uipo.^ 

\ 

XY 

E(Xj-X/)(Vk-VkF) = 


(4.4) 
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where tie oj^'s are known, and let I be the (N^ax + f\nax) * (Nmax + %tx) 
covariance matrix with entries Let 2 the me|suremegt 

errors of (u^, v<j, be independent with variances wgOg , u^Oy , (itgo^ . 
Suppose all random variables are normally distributed. Then 

Theorem; (6. Wahba, in preparation) The Bayes estimate 

A A A 

® ® (X : Y) of 6 a (X : Y) is the minitnizer of 


1 

n 


n 

I 

i=l 


Ui 


h - 

J ^ 


^ax f. 

I »kH|tS(p,) 

k=l 


1 

a 


(4.5) 


t «>((x-xf; v-yf) I-> (x-xf; y-yf)') 



where 



0 ) a uip/upn, X = u>o/bjn, 6 = bj/b^ . 


Since the expression to be minimized i: a quadratic form in the components 
of (X:Y), (for given id, X, 6), the minimizer can be readily expressed as 

the solution to a linear system. It is conjectured that the bandwidth 
parameter X, the partitioning parameter 6 and the "error balancing" parameter 
ID, can all be estimated (simultaneously!) from the data by GCV. Note that a 
large id suppresses forecast relative to observational data and a large 6 
suppresses fast modes relative to slow modes. The results reported in 
Section 3 suggest that at least one bandwidth and one (appropriately chosen) 
partitioning parameter can be chosen by GCV. Part of the art of this 
approach is to choose these "tuning" parameters so that they are (a) the 
important ones and (b) their determination is relatively "well posed." 
Sensitivity of optimum interpolation to bandwidth parameters has recently 
been discussed by Hollingsworth (1982), Lorenc (1981) and others. 

To be practical it may be necessary to assume an oversimplified structure 
for T and/or to break up the variational problem into a series of sub- 

A 

problems. The usual Kalman filtering theory would give 6 as the minimizer of 
the above expression, with X - 0 and with (wpl), the forecast error 
covariance, given by a recursion formula in time. Our description above. 
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Impllcltly assumes the existence of a limiting (upi:). The forecast of ^11 
et al. would, roughly speaking reduce In this context , to setting X ■ 0 and 

A 

constraining Y to be 0. The above theorem generalizes to the simultaneous 

analysis of all levels. Then satellite radiance data (which "cuts across" 
all levels) can be Included as another term In the variational formulation. 
Physical constraints can. In principle be Included as side conditions. It 
also appears that non-linear balancing constraints related to those discussed 
In, e.g. Tribbla (1982), whose purpose Is to minimize the propagation of 
unwanted gravity waves In the non-linear forecast equations, may also be 
Incorporated as either weak or strong constraints. 
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